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A MODIFIED GALERKIN / EINITE ELEMENT METHOD FOR 
THE NUMERICAL SOLUTION OE THE 
SERRE-GREEN-NAGHDI SYSTEM 

DIMITRIOS MITSOTAKIS, COSTAS SYNOLAKIS, AND MARK MCGUINNESS 


Abstract. A new modified Galerkin / Finite Element Method is proposed 
for the numerical solution of the fully nonlinear shallow water wave equa¬ 
tions. The new numerical method allows the use of low-order Lagrange finite 
element spaces, despite the fact that the system contains third order spatial 
partial derivatives for the depth averaged velocity of the fluid. After studying 
the efficacy and the conservation properties of the new numerical method, we 
proceed with the validation of the new numerical model and boundary condi¬ 
tions by comparing the numerical solutions with laboratory experiments and 
with available theoretical asymptotic results. 


1. Introduction 


The motion of an ideal (inviscid, irrotational) fluid bounded above by a free 
surface and below by an impermeable bottom is governed by the full Euler equations 
of water wave theory, [63]. Because of the complexity of the Euler equations, a 
number of simplified models describing inviscid fluid flow have been derived such 
as various Boussinesq type (BT) models. The Serre-Green-Naghdi (SGN) system 
can be considered to be a BT model that approximates the Euler equations, and 
models one-dimensional, two-way propagation of long waves, without any restrictive 
conditions on the wave height. The SGN system is a fully non-linear system of the 
form, 

ht + {hu)x = 0 , (la) 

[h + 7),^] ut + gh{h -j- b)x + huu^ + Q!^u + Q^u = 0 , (lb) 

where 


h{x,t) = r]{x,t) — b{x) , 


(Ic) 


is the total depth of the water between the bottom b(x) and the free surface elevation 
i](x,t), u(x,t) is the depth averaged horizontal velocity of the fluid, and g the 
acceleration due to gravity. The operators and Q!^ depend on h, and are 

defined as follows: 


Th^w = h 


hx^x “f ^hbxx 


W-- [h^Wx]^ , 


(Id) 




1 

3 


[h^iwwxx -wl)]^ , 


(le) 
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Qb'w = 2 + wWxbx)]^--^h^ [wwxx - wl] bx+hup'bxbxx+hwwxbl . 

(If) 


In dimensional and unsealed form, the independent variable, x G R is a spatial 
variable and t > 0 represents the time. 

The SGN equations as derived by Seabra-Santos et.al. in have also been 
derived in a three-dimensional form in [36] and in a different formulation by Green 
and Naghdi m- In the case of a flat bottom (i.e. = 0) ([1]) is simplified to the 

so-called Serre system of equations derived first by Serre EIllSl] and re-derived 
later by Su and Gardner, [5S]. For these reasons the equations m are also known 
as the Serre, or Green-Nagdhi, or Su and Gardner equations. We will henceforth 
refer to them here as the Serre-Green-Naghdi (SGN) equations. 

Under the additional assumption of small amplitude waves (i.e. the solutions 
are of small amplitude), the SGN system reduces to the Peregrine system, [T7] : 


ht + {hu)x = 0 , 

(2a) 

h 

9\^ “1” “1“ UUx 2 g '^xxt — 0 ■ 

(2b) 


Peregrine’s system belongs to the weakly dispersive and weakly nonlinear BT sys¬ 
tems. There are also other BT systems that are asymptotically equivalent to Pere¬ 
grine’s system, cf. [431144] . Differences between the SGN equations and BT models 
are explained in |19) . In the same work the inclusion of surface tension effects have 
been included and explained in detail. 

Although equations ([U can be derived from the SGN system, their solutions have 
different properties. For example the solutions of the SGN equations are invariant 
under the Galilean boost, while the respective solutions of ([2]) are not. It is noted 
that the SGN equations have a Hamiltonian formulation, [39l [31] . Specifically, for a 
stationary bathymetry the SGN system conserves the total energy functional m- 


I{t) = / gr]"^ + hu^ + T^u ■ u dx , (3) 

Jr 

in the sense that I{t) = 7(0), for all t > 0. The conservation of this Hamiltonian will 
be used to measure accuracy and conservation properties of the proposed numerical 
methods. 

Although both systems are known to admit solitary wave solutions propagating 
without change in their shape over a horizontal bottom y = —bg, only the solitary 
waves of the SGN system have known formulas in a closed form. Specifically, a 
solitary wave of the SGN system with amplitude A can be written in the form: 


hs{x,t) = bo + Asech^ [A(x — Cgt)], Us (x, t) = cq ( 1 — 


h{x, t)) 


A = 


3A 


4&g(6o + A) 


. I ^ 

and Cs = co\/ 1 + — , 


(4a) 

(4b) 


where Cg is the phase speed of the solitary wave and cg = i/gfog is the linear wave 
speed. 

In addition to the above properties the system o is known to have several 
favourable well-posedness properties. For example, it is locally well-posed in time 
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for smooth bottom functions in R. Specifically, if C^{I) is the space of bounded and 
continuously differentiable functions on the interval I and (/) denotes the usual 
Sobolev space of s-order weakly differentiable functions on I, and if 6 £ C'^(R), 
K > 1/2, s > K + 1 with the initial condition is (?7o,'ao) € x then 

there is at least a maximal time Tmax > 0, such that the SGN equations admit 
a unique solution (??, u) € H^(R) x cf. [31]. We also mention that the 

solution will satisfy a non-vanishing depth condition when the initial data satisfy 
the same non-vanishing condition, i.e. there is an a > 0 such that for alH > 0 

y(x, t) — b{x) > a > 0 . (5) 

For more information on the model equations, including the derivation, theory and 
justification, refer to [551IM1 I5]. 

We consider here the SGN equations on a finite interval I = [a, 6]. In addition to 
periodic boundary conditions, we will study physically relevant reflective boundary 
conditions. The linearized equations © about the trivial solution coincide with 
the linearized Peregrine system. Initial-boundary value problems for Peregrine’s 
system have been studied in [26] . where the existence of solutions was proven for 
boundary data for u{a, t) and m(6, t). In a similar manner, the boundary conditions 
u{a, t) = u{b, t) = 0 for t > 0 are sufficient for the system (|T]) to model the reflection 
of the waves on solid wall boundaries. These specific wall boundary conditions 
have been used widely to describe reflection of waves on the boundaries for various 
numerical models including BT systems such as the Nwogu system, [34] . 

The numerical discretization of the SGN equations is a challenging problem 
due to their complicated form. Not only do the operators in front of the temporal 
derivatives depend on the unknown function h(x, <), but they also contain high order 
derivatives in nonlinear terms. Because of the presence of those high order spatial 
partial derivatives the solution to the system should be smooth. Recently several 
schemes have been proposed such as Finite Difference / Finite Volume Schemes 
and Discontinuous Galerkin methods [351 HSUS). Although 
these methods are very useful to study practical problems, such as the runup of 
nonlinear waves on slopes, they can be highly dissipative methods in the sense that 
they introduce numerical dissipation or dispersion, due to the approximation of 
the nonlinear terms by dissipative flux functions and the use of low order slope 
limiters. These methods appear to have good conservation properties though when 
high-order WENO methods are used or other non-classical numerical fluxes are 
used, [37] . 

Some other highly accurate numerical methods have been developed for the Serre 
equations in the case of a horizontal bottom, such as spectral methods [23] and 
Galerkin / Finite Element Method (FEM) [35]. Although these methods appear to 
have satisfactory conservation properties, it is very difficult to extend them to the 
SGN equations especially in the case of two horizontal dimensions. Eor example 
the standard Galerkin method of [42] requires tensor products of cubic splines in 
order to be consistent in two spatial dimensions in a similar manner to |22j . 

In this study, a new modified Galerkin / Finite Element Method (EEM) is pro¬ 
posed for the numerical solution of the SGN equations. This method allows the 
use of low-order finite elements such as piecewise quadratic (P^) or even piecewise 
linear (P^) Lagrange finite elements. Two of the main advantages of this method is 
that it is highly accurate, and it has very good conservation properties. Other ad¬ 
vantages are the sparsity of the resulting linear systems, the low complexity of the 


4 DIMITRIOS MITSOTAKIS, COSTAS SYNOLAKIS, AND MARK MCGUINNESS 

algorithm due to the use of low-order finite element spaces and finally its potential 
to be extended to the two-dimensional model equations of [35]. 

Similar techniques that reduce the requirement of high-order finite elements (for 
example the use of cubic splines) have been used previously for weakly-nonlinear 
Boussinesq systems in |611I621[2T| where the second derivative in the linear dispersive 
terms has been replaced by either the discrete Laplacian operator or the solution of 
an intermediate problem. In the case of the Bona-Smith type of Boussinesq systems 
with wall boundary conditions, the modified Galerkin method converges at an op¬ 
timal rate showing great performance contrary to the suboptimal convergence rates 
achieved with the standard Galerkin / FEM method, [21]. Similarly to the behav¬ 
ior of the modified Galerkin method for the weakly-nonlinear Boussinesq equations, 
the proposed FEM scheme for the SGN equations can achieve optimal convergence 
rates depending on the choice of the trial function spaces, contrary to the subop¬ 
timal rates obtained for the standard Galerkin method for the SGN system and 
also for Peregrine’s system as shown in [J. Specifically, in order to achieve optimal 
convergence properties one may use spaces of piecewise linear elements for the free 
surface elevation and piecewise quadratic elements for the horizontal velocity. 

The validated numerical method is applied to study the SGN equations with 
reflective boundary conditions in a systematic way through a series of numerical 
experiments. In particular, we focus on the following issues: 

• accuracy of the modified Galerkin method and invariant conservation; 

• reflection of solitary waves at a vertical wall; and 

• shoaling of solitary waves on plain or composite beaches. 

The convergence properties of the new numerical method are also tested in the case 
of periodic boundary conditions. For more information about the behavior and the 
properties of the Galerkin / FEM method with periodic boundary conditions we 
refer to [42]. 

This paper is organized as follows. Section [5] presents the fully discrete schemes 
for the SGN equations. In Section [3] we study the convergence, accuracy and 
numerical stability of the modified Galerkin method. Finally, Section |4] presents 
computational studies of shoaling and reflected waves validating both the choice of 
boundary conditions and the numerical scheme. We close the paper with conclu¬ 
sions in Section |5j 


2. The numerical methods 

We consider the initial-boundary value problem (IBVP) comprising System ([T]), 


subject to reflective boundary conditions: 

ht + {hu)x = 0 , (6a) 

\h T^'\ ut + gh{h b)x + huux + Q!^u + Q!f,u = Q , (6b) 

u{a, t) = u(b, t) = 0, h(x, 0) = ho(x), u{x, 0) = uo{x) , (6c) 

where again 

= h hxbx + ^hbxx + bl [h^Wx] ^ , (6d) 

Q’^w = -i [h^iwwxx -Wx)]^ , (6e) 




A MODIFIED GALERKIN METHOD FOR THE SGN SYSTEM 


5 


QbW = ^ [h'^iw'^bxx + wwxbx)]^-^h‘^ [wwxx - wl] bx+hw'^bxbxx+hwwxbl , 

(6f) 

X £ {a, b) C M. and t G [0,T]. We assume that ([6]) possesses a unique solution, 
such that h and u are sufficiently smooth and, for any t G [0, T], in suitable Sobolev 
spaces 

u{x,-)£H^+\ 

where s > 1. Here and below, || • ||s denotes the standard norm in while Hq will 
denote the subspace of iJ® whose elements vanish at a: = a and x = b. We also use 
the inner product in denoted by (•,•), which is 

fb 

{u,v) = uv dx , 

J a 

for all u,v £ L^. 

A spatial grid of the interval [o, b] is a collection of points Xi = a + i Ax, for 
i = 0,1, ■ ■ ■ ,N, where Ax is the grid size, and N £ N, such that Aa: = {b — a) /N. 
Let {h, u) £ Sh X Su be the corresponding spatially discretized solutions of the 
Galerkin / Finite Element Method (FEM) for suitable finite-dimensional spaces Sh 
and Su- First we present the standard Galerkin / FEM semidiscretization. 


2.1. The standard Galerkin method. For the standard Galerkin method, we 
consider the space of smooth splines 

= {^£C^-^[a,b]\ ^1 G P^ 0<z< TV-l} , 

where P'’ is the space of polynomials of degree r. The standard Galerkin method 
requires r > 3. Here, we take r = 3. The trial function space for the first equation 
and for the solution h is chosen as Sh = 5'’’, while for the second equation and 
for the approximation of the depth averaged velocity of the fluid u is Su = S"^ C\ 
G (^[a, &]| (p{a) = ip{b) = 0}. To state the associated semi-discrete problem, let 
4> £ Sh and Ip £ Su he arbitrary test functions. After taking inner products, and 
using integration by parts, the semi-discrete problem takes the form: 


{ht,4') + 


,</> =0 , 


B{ut,ip;h) + (ji g{h + b)x+uu,, 


(7a) 

+ Q(u, -0; h) + Qb{u, Ip-,h) = 0 , 

(7b) 


where B, Q and Qb are defined for any u!,ip £ Su as 

+ ^(h^Wx,'ipx'j , (7c) 

Q{u}, tp-, h) {h^ [wujxx - i^l] , , (7d) 

Qb(^^, 0? h) = — (cU bxx “t“ bJUJxbx') ^ 

hbx |h [ujbJxx - - bJ%xx - ujuJxbx'j , 0^ . (7e) 

This is a system of ordinary differential equations (ODEs). Given the initial 
conditions h{x,0) = ho = Vh{ho) and M(a:,0) = uq = Vu{uo) where Vh and Vu are 



B{uj, ip-,h) = i h 


1 -f hxbx -f ^hbxx 


a ;,0 
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appropriate projections on Sh and Su respectively, we assume that the system © 
has a unique solution. Appropriate projections of the initial conditions could be 
the standard L^-projections on Sh and Su, defined as Hq G Sh and Uq € Su such 
that 


Hod) dx = 


ho4> dx, for all (j) & Sh 


and 


IIqi/) dx= uoi) dx, for all ip G Su 


The presence of the term lOxx in Q and Qb requires the use of at least (7^ smooth 
splines. In particular, one may use cubic splines, which correspond to Sh = S^, 
i.e., r = 3. A basis for the space Sh for a uniform grid of mesh-length Ax can be 

j = —1,0, • • • , -I- 1 where 


[a,b] 


formed by the functions = B[x — Xj/Ax) 

X-i = a — Ax, xn+i = b Ax and 

K^ + 2)^ 

t[ 1 -l- 3(a; -|- 1) -|- 3(a; -I- 1)^ — 3(a; -|- 1)^] 

in I I Qi'i Q/1 .„\3i 


B{x) = < 


^[1 + 3(1 - a;) -h 3(1 - a;)^ - 3(1 - a;)^] 
*(2-x)3 


-2 < X < -1 
-1 < X < 0 

0 < X < 1 
1 < X < 2 
X GR-[-2,2] 


The basis of Su = Sq can be described by the functions ipj (x) = (pj (x) for 2 < j < 
N — 2, plus four functions 'tpQ,ipi,ipN-i,'tpN, taken as linear combinations of the 
(p-i,(po,(pi and (pN-i,(pN,(pN+i, which are such that ipoia) = ipi{a) = ipN-i{b) — 
i^N{b) = 0. For example, we take ipo = 4’o — ^(p-i and ipi = (pi — (p-i, [IH]- The 
convergence properties of the standard Galerkin method with cubic splines are very 
similar to those of the modified Galerkin method with S^ elements and thus are 
not presented here. Some details and numerical experiments with the standard 
Galerkin method can be found in |32] . 


2.2. The modified Galerkin method. In real-world applications, the use of low- 
order finite element methods can lead to faster computations. We derive a numerical 
method that does not require high-order finite element spaces. This can be done 
by using similar techniques to those proposed in ismsmi], but applied to the 
nonlinear term uuxx- The resulting method is a modified Galerkin method that 
allows the use of Lagrange elements as low as in order, consisting of piecewise 
linear functions. We proceed with the derivation of the modified Galerkin method 
but first we introduce the notation for the Lagrange finite element spaces that are 
subspaces of It is noted that we cannot use Lagrange finite elements with the 
standard Galerkin method due to the second order derivative. 

The space P’’ of Lagrange finite elements is defined on the grid x^ = xq -I- iAx, 
i = 0,1,-- ■ ,N as: 

P^ = {x&C°[a,b]\x\[x,,x,^,]GF\ 0<i<N-l} , (8) 

We will restrict the analysis below to P^, P^ and P^ Lagrange finite element spaces. 
The P^ Lagrange elements can be defined using the basis functions 

f if X G [xi_i,Xi] , 

Xi{x) = < if X G [xi,Xi+i] , 

[ 0 , otherwise . 
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The Lagrange finite element space can be defined on the same grid Xi = xo+i^x 
for z = 0,1, • • • ,N and at the midpoints aii+ 1/2 = Xi + Ax/2 for z = 0,1, • • • , iV — 1 
using the basis functions 


and 


with 


and 


X^{x) = 
Xz+l/2{x) ='4’ 


X — Xi 


Ax 
X — \ 

Ax J 


, z = 0,1,-- - , 


, i = ,iV-l 


(1 + x)(l + 2x), —1 < a; < 0 , 
ip{x) = ■^ (1 — x)(l — 2x), 0 < X < 1 , 

0 otherwise , 


il){x) = 


1 — 4x^, |x| < 1/2 , 

0, otherwise . 

Then a function w € is written as 

N N-1 

w{x) = E w{x^)Xi{x) + ^ w{x^+il2)Xt+l/2{x) . 


i=0 


i=0 


For the construction of the Lagrange basis function of the general space P'' we refer 
to [24]. We also consider the Lagrange finite element spaces Pq = {x G P’'|x(a) = 
x{b) = 0}. These spaces are subspaces of iLg and will be used to approximate 
the depth averaged horizontal velocity of the water. For example, we will use the 
spaces Sfi = P^ and Su = Pq, for some integers r,q or we will use the spaces of 
smooth splines described in the previous section. For more information related to 
Lagrange finite element spaces and its approximation properties, we refer to [24] . 
It is worth mentioning that, for the discretization of the bottom boundary, we use 
the projection of the bathymetry and its derivatives. 

Given v € PIq, we define the non-linear discrete Laplacian operator : PIq —^ Su 
such that 


{d'^v, -0j = - {vl, 0 ) - {vvx,i/x), for all ip € Su ■ 


(9) 


In fact, the function d^v G Su approximates the function vvxx as if v was a smooth 
function and substitution in (l7d]) and (ITel) leads to the modified Galerkin semi¬ 
discretization: 


{ht,(j))i^(hu)x,4>j = 0, (p G Sh , (lOa) 

B{ut,'ip;h)(h gih-Gb),^-Guux , 0 ^ + 2(11, 0 ; h) + 0 ; iz) = 0, ip G Su 

(10b) 


where B, Q and Qb are defined for any w, 0 S 5'^ as 


B{u!, ip]h) = I h 


1 P hxbx + —hbxx + bi 


, 0 ) + ^ {h^Wx,ipx') 


Q{uj,tp;h) = i (h^ 


d^LO — I 


,0x) 


(10c) 

(lOd) 
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- 2 bxx ) ^ ‘4^X^ 


i {hhx [^' 


— ut 


^ ^xx ^^x^x 


(lOe) 


Although we assume that the bottom function has the appropriate smoothness, 
the absence of second order spatial derivatives of the depth integrated horizontal 
velocity in the semidiscrete scheme allows the use of Lagrange elements, as well 
as high-order elements such as cubic or quintic splines. In what follows, we test 
the efficiency of the modified Galerkin method using the spaces P^, and of 
Lagrange elements and the space of cubic splines with periodic and reflective 
boundary conditions. 


Remark 2.1 (Mass lumping). In order to compute the non-linear discrete Lapla- 
cian dw one needs to solve the linear system obtained by the discretization of m- 
In the case of wall boundary conditions if, for example, ifi denotes basis functions 
of Su the system can be written as A4w = f where the mass matrix is a banded 
matrix with entries Aiij = (pjjijipj) and fi = — ~ i'i’Vx,'if>j)- To improve 

the speed of the numerical method, one may apply the method of mass lumping in 
the formulation of the matrix M.. This can be done, for example in the case of 
quadratic Lagrange elements, by approximating the integrals with Simpson’s rule. 
This leads to a diagonal matrix that can be inverted trivially. All the numerical 
experiments with P^ and P^ elements have been performed with mass lumping, in 
addition to the standard matrix formulation with comparable results. 


Remark 2.2. The choice of the discrete Laplacian is not unique. For example, 
when periodic boundary conditions are used, then one may consider the linear dis¬ 
crete Laplacian, which replaces the term Uxx o,s proposed in [smsiiii]. 

Remark 2.3. The presence of the term bxx in the semidiscrete scheme implies the 
typical requirement of a smooth bottom. When a piecewise linear bottom topography 
is given, then the use of an appropriate projection onto the finite element space 
(such as the elliptic projection) or local smoothing of the bottom is required. In 
this paper the bathymetry is usually a piecewise linear function and thus we use the 
smoothing method described in [1]. 

Remark 2.4. The modified Galerkin method can be expressed in an equivalent 
way by introducing an additional independent variable w satisfying the equation 
W = UUxxt EH- This increase of the degrees of freedom in the model equations by 
one allows the approximation of the unknown function u by low-order finite element 
spaces such as P^ or P^ finite element spaces. 

2.3. Temporal discretization. In [42], it was shown by numerical means that 
the standard Galerkin method for the discretization of the SGN equations with flat 
bottom leads to a system of ODEs that is not stiff. Also, the classical, explicit, four- 
stage, fourth order Runge-Kutta (RK) method, described by the following Butcher 
tableau: 


0 

0 

0 

0 

1/6 

1/2 

0 

0 

0 

1/3 

0 

1/2 

0 

0 

1/3 

0 

0 

1/2 

0 

1/6 

0 

1/2 

1/2 

1 



(11) 
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is able to integrate the respective semi-discrete system numerically in time with¬ 
out imposing restrictive stability conditions on the ratio At/Ax, but only mild 
restrictions on the mesh length such as At/Ax < 2 for smooth solutions. 

Concerning the time integration, the modified Galerkin method has very similar 
behavior to the standard Galerkin method. Upon choosing appropriate basis func¬ 
tions for the spaces Sh and Su, the semidiscrete system m represents a system of 
ODEs. We use a uniform time-step At such that At = T/K ior K G N. The tem¬ 
poral grid is then = nAt, where n = 0,1, • • • ,K. Given the ODE y' = $(t, y), 
one step of this four-stage RK scheme (with y" approximating y{t")) is: 
for * = 1 —4 do 

f = y^ + AtY:~=\a,,y^’^ 
yn,t _ evaluated at t"’* = t^+ nAt 

end for 

where Oy , n, h are given in table (fTT|l . Applying this scheme to (fT0)l and denoting 

by iJ" and C/" the fully discrete approximations in Sh and Su of 

respectively, leads to Algorithm [TJ Given bases {c/i} of Sh and {tpi} of Su, the 


Algorithm 1 Time-marching FEM scheme for the IBVP of the system (nni) 
= r{ho} 
u° = V{uo} 
for n = 0— Ido 
for i = 1 ^ 4: do 


W = U" -b AtY/yJiarjU'^'^ 

V') = evaluated at t”’® = -b nAt 




(h^ \g{W - b), + WUl 


P{U\(/-H\b)=0 


where V{U\ (/; H\ b) = Q{U\ (/; W) + Qi{U\ </>; Hp 

end for 


jjn+l = 

end for 


implementation of Algorithm [T] requires solving at each time step the following 
linear systems: 

(a) Four linear systems with the time-independent matrix {(/i, (/j); 

(b) Four linear systems with the time-dependent matrix /i); 
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(c) Four linear systems with the time-independent matrix . 

The four linear systems in (a) arise from the discretization of the equation (llOaL 
while the four linear systems in (c) arise from the computation of the discrete 
Laplacian All of these matrices are banded and symmetric, and one can use 
either direct methods for banded systems or classical iterative methods for sparse 
systems. To approximate the inner products, we use Gauss-Legendre quadrature 
with 3 nodes per Ax for elements, 5 nodes for and P^, and 8 nodes for 
elements. It is noted that most of the experiments with P^, P^ and mixed 
elements have also been tested using adaptive time-stepping methods such as the 
Runge-Kutta-Fehlberg method, [5^, to ensure that the errors introduced by the 
temporal integration are negligible. 

3. Accuracy and convergence 

As it was pointed out in the introduction, the ‘classical’ Boussinesq (cB) system 
is similar in structure to the SGN system, given that both systems admit the same 
number of boundary conditions, and the equation for the free surface is the same. 
We thus expect similar behavior for the convergence of the numerical method. In 
[2], it was shown that the convergence of the Galerkin / FEM method for the cB 
system with periodic boundary conditions is optimal, while in nig it was shown 
that the convergence of the Galerkin / FEM method with cubic splines for the cB 
system is suboptimal, and that in order to achieve an optimal rate of convergence, 
a non-standard method should be used. Specifically, it was shown that for wall 
boundary conditions in the case of P^ elements, the following error estimates hold: 

max \\h— h\\ < C max llu — u|| < C Ax^ , (12) 

0<t<T 0<t<T 

which are suboptimal for h and optimal for u. In the case of cubic splines, the 
following error estimates hold: 

max ||/i — hll < C Ax^'^-i/ln ——, max llu — ull < C Ax^\/ln —— , (13) 

o<<t" V Ax o<t<T" V Ax ^ ^ 

which is suboptimal in both h and u but the factor In Ax is not dominant and 
generally it has an effect on the accuracy of the method that is too small to observe. 
In [32) . it is shown, by numerical means that the standard Galerkin method with 
cubic splines for the SGN system with wall boundary conditions appears to have 
similar convergence properties. 

The standard Galerkin method for the initial-periodic boundary value problem 
for the SGN system has been studied in [42], where it was shown that the con¬ 
vergence is optimal. It is worth mentioning that the approximation m can be 
used also with periodic boundary conditions, and the resulting modified Galerkin 
method has optimal convergence rates including all the advantages of the modifi¬ 
cation. Specifically, the errors in the norm for the periodic case are of O(Ax^) 
for the P^ elements for both h and u, O(Ax^) for the P^ elements and for the 
elements is of 0(Ax'^). Because the results for the periodic problem are very 
similar to those of [35] and to the modified Galerkin method for the SGN system 
with periodic boundary conditions, we don’t present them here. 

We continue with the wall boundary conditions. In order to study the conver¬ 
gence of the numerical method, we use the SGN equations in nondimensional but 
unsealed form. We start with the IB VP ([6|) . Because we don’t know any analytical 
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solutions for this problem, we consider the non-homogenous problem, which has 
the exact solution h{x,t) = 1 + e^*(cos(7ra:) + x -\-2) and u{x,t) = x sin{TTx) 
satisfying the equations (I5al) and (l6bl) with the appropriate right-hand sides. The 
solution is computed using the modified Galerkin method in the interval [0,1] and 
for t G (0,T] with T = 1. The finite element spaces of our choice are and 
Lagrange finite element spaces and the space of cubic splines. It is noted that 
we also tried other “exact” solutions with different right-hand sides in (I6a1l and 
(HEl) to verify the computed rates of convergence. The results have been all very 
similar, and thus we only present the results of one case. 

In order to study the accuracy and the convergence of the modified Galerkin 
method, several error indicators have been computed. The computed errors are 
normalized and defined as 

\\F{x,T;Ax) T)\\s 


Es[F] = 


(14) 


ll-^exact(2^; "^)lls 

where F = F{-;Ax) is the computed solution, i.e., either F[ k. h{x,T) or ~ 
u(x, T), inexact corresponding exact solution and s = 0,1, 2, oo correspond to 

the F[^, FI^ and L°° norms, respectively. The analogous rates of convergence 
are defined as 


rate for Es[F] 


HE,[F{-; Axk-i)]/E,[F{-; Axfc)]) 

ln(Aa;fc_i/Aa;fe) 


(15) 


where Axk is the grid size listed in row k in each table. To ensure that the errors 
incurred by the temporal integration do not affect the rates of convergence we use 
At ^ Ax while we take Ax = 1/N. 

Tables (U [2] and [3] present the errors and the corresponding rates of convergence 
of the modihed Galerkin method with P^ finite elements. It is shown that the 
rate of convergence is suboptimal for the total depth h and optimal for velocity u. 
Specihcally, Table [T] suggests that \\h — h\\ — 0{Ax^^‘^) and ||m — i2|| = O(Ax^). 
Table [2] suggests that ||/i —h||i = 0(Ax^^^) and ||m —u||i = 0(Ax). Finally, the L°“ 
estimates shown in[3]are worse than the estimates, suggesting ||/i—h|loo = 0(Ax) 
and ||m — u||oo = O(Ax^). All these results are very similar to those obtained 
theoretically and numerically in the case of Peregrine’s system [4]. 


N 

Eo[ff] 

rate for Eq[F[] 

Eo[U] 

rate for Eq\U] 

10 

1.4661 X 10-^ 

- 

2.9141 X 10-^ 

- 

20 

3.3761 X 10-3 

2.1186 

2.3778 X 10-3 

3.6154 

40 

1.1967 X 10-3 

1.4963 

5.6477 X 10-4 

2.0739 

80 

4.4536 X 10-* 

1.4260 

1.3856 X 10-4 

2.0271 

160 

1.6179 X 10-4 

1.4609 

3.4262 X 10-3 

2.0159 

320 

5.7982 X 10-3 

1.4804 

8.5165 X 10-3 

2.0083 

640 

2.0638 X 10-3 

1.4903 

2.1229 X 10-3 

2.0042 


Table 1. Spatial errors and rates of convergence for the exact 
solution with P^ finite elements using the norm. 


In the case of quadratic P^ elements, the results with the same initial conditions 
are similar, but the convergence rates for h are different to the expected suboptimal 
rates, analogous to the P^ case. The respective convergence rates for u are optimal 
as expected. Tables HI [5] and |6] present the errors and the convergence rates for the 
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N 

Ei[H] 

rate for Ei[P[] 

Ei[U] 

rate for Ei \U] 

10 

2.8836 X 10-1 

- 

1.6945 X 10-1 

- 

20 

1.8855 X 10-1 

0.6129 

6.6788 X 10-2 

1.3433 

40 

1.3684 X 10-1 

0.4625 

3.2009 X 10-2 

1.0611 

80 

1.0135 X 10-1 

0.4331 

1.5859 X 10-2 

1.0131 

160 

7.3388 X 10-2 

0.4658 

7.9074 X 10-3 

1.0041 

320 

5.2503 X 10-2 

0.4831 

3.9492 X 10-3 

1.0016 

640 

3.7340 X 10-2 

0.4917 

1.9736 X 10-3 

1.0007 


Table 2. Spatial errors and rates of convergence for the exact 
solution with finite elements using the norm. 


N 

E^[H] 

rate for Eoc:[H] 

Eo.[U] 

rate for Eoc\U] 

10 

3.4775 X 10-2 

1.4587 

3.9193 X 10-2 

1.4068 

20 

1.1039 X 10-2 

1.6554 

5.7364 X 10-3 

2.7724 

40 

5.8666 X 10-3 

0.9121 

1.0352 X 10-3 

2.4701 

80 

3.4929 X 10-3 

0.7481 

2.1191 X 10-4 

2.2884 

160 

1.8861 X 10-3 

0.8890 

5.2158 X 10-3 

2.0225 

320 

9.7698 X 10-4 

0.9490 

1.3480 X 10-3 

1.9520 

640 

4.9684 X 10-4 

0.9755 

3.4245 X 10-3 

1.9769 


Table 3. Spatial errors and rates of convergence for the exact 
solution with finite elements using the norm. 


modified Galerkin method with P^ elements. Table S] suggests that the \\h — h\\ = 
O(Ax^) and ||u — u\\ = 0{Ax^). Table [5] suggests that \\h — ii||i = 0{Ax) and 
||u — fill = 0{Ax^). Finally, Table |6] suggests the same estimates in the L°° norm 
with norm, i.e. \\h — ^||oo = O(Ax^) and ||u — u||oo = O(Ax^). The same 
behaviour was observed when we used P^ elements, i.e. the convergence rate for 
h was suboptimal \\h — h\\s = 0{Ax^~^) while the convergence for u was optimal 


U\\s = 

0(Aa;4 ®), for s 

= 0,1. 



N 

Eo[H] 

rate for Eo[H] 

Eo[U] 

rate for Eo[U] 

10 

3.1489 X 10-3 

- 

1.5559 X 10-3 

- 

20 

6.1429 X 10-4 

2.3579 

1.1718 X 10-4 

3.7310 

40 

1.3897 X 10-4 

2.1441 

9.8909 X 10-3 

3.5665 

80 

3.3909 X 10-3 

2.0351 

1.0744 X 10-3 

3.2025 

160 

8.4285 X 10-3 

2.0083 

1.2858 X 10-^ 

3.0628 

320 

2.1019 X 10-3 

2.0035 

1.5874 X 10-8 

3.0180 

640 

5.2473 X 10-^ 

2.0021 

1.9773 X 10-® 

3.0051 


Table 4. Spatial errors and rates of convergence for the exact 
solution with P^ finite elements using the norm. 


We also studied the errors and convergence rates of a mixed modified Galerkin 
method, where we considered St = P^ and S'„ = P^ i.e. we used linear elements 
for the approximation of h and quadratic elements for the approximation of u. 
This was the only case where optimal rates of convergence for both h and u were 
obtained. Tables [THSl suggest that ||/i—ii||s = and ||u —{t||s = 0{Ax^~‘^) 
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N 

Ei[H] 

rate for Ei[H] 

Ei[U] 

rate for Ei [U] 

10 

1.6038 X 10-^ 

- 

1.8439 X 10-2 

- 

20 

7.6716 X 10-2 

1.0640 

3.3234 X 10-3 

2.4721 

40 

3.7406 X 10-2 

1.0362 

6.4326 X 10-* 

2.3692 

80 

1.8545 X 10-2 

1.0122 

1.4083 X 10-* 

2.1915 

160 

9.2498 X 10-3 

1.0036 

3.3536 X 10-3 

2.0701 

320 

4.6179 X 10-3 

1.0022 

8.2655 X 10-® 

2.0206 

640 

2.3064 X 10-3 

1.0016 

2.0585 X 10-® 

2.0055 

Table 5. Spatial errors and rates of convergence 
solution with P^ finite elements using the norm. 

for the exact 


N 

Eoo[H] 

rate for Eac[H] 

Eo.[U] 

rate for Eoc [U] 

10 

7.4168 X 10-3 

- 

2.7712 X 10-3 

- 

20 

1.7953 X 10-3 

2.0466 

2.4407 X 10-4 

3.5051 

40 

4.0253 X 10-'* 

2.1571 

2.3644 X 10-3 

3.3678 

80 

1.0353 X lO-'* 

1.9589 

1.8056 X 10-3 

3.7109 

160 

2.5922 X 10-3 

1.9979 

1.6493 X 10-^ 

3.4525 

320 

6.5561 X 10-3 

1.9833 

2.0640 X 10-8 

2.9984 

640 

1.6390 X 10-3 

2.0000 

2.5810 X 10-3 

2.9994 


Table 6. Spatial errors and rates of convergence for the exact 
solution with finite elements using the L°° norm. 


for s = 0,1 while in the maximum norm we found that ||/i — /i||oo = 0{Ax^) 
and ||u — m||oo = 0{Ax^). Using mixed elements of higher order does not give 
analogously optimal results. For example, when we took Sh = and = P^ 
the convergence was suboptimal for h and optimal for u. More specifically, the 
errors were \\h — h\\s = 0(Aa;^“®) and ||m — m||s = 0(Aa;^“®) for s = 0,1 while in 
the maximum norm we found that \\h — h\\oa = 0{Ax^) and ||m — m||oo = 0{Ax*). 
Observe that the convergence for h achieved with the mixed P^ — P^ elements was 
the same with the mixed P^ — P^ elements. 


N 

Eo[H] 

rate for Eo[H] 

Eo[U] 

rate for Eo[U] 

10 

1.3308 X 10-3 

- 

1.4064 X 10-3 

- 

20 

2.6812 X 10-4 

2.3113 

9.9809 X 10-3 

3.8167 

40 

6.2091 X 10-3 

2.1104 

9.1910 X 10-3 

3.4409 

80 

1.5350 X 10-3 

2.0161 

1.0541 X 10-3 

3.1241 

160 

3.8219 X 10-3 

2.0059 

1.2825 X 10-'^ 

3.0390 

320 

9.5358 X 10-'^ 

2.0029 

1.6047 X 10-8 

2.9986 

640 

2.38157 X 10-^ 

2.0014 

1.9751 X 10-3 

3.0223 


Table 7. Spatial errors and rates of convergence for the exact 
solution with Sh = P^ and Su = P^ using the norm. 


Finally, Tables [TOl E] and M present the errors and convergence rates for the 
modified Galerkin method with cubic splines, i.e. with elements. In this case, the 
results are similar to those of the standard Galerkin method for the SGN and Pere¬ 
grine’s system. Specifically, Table [TUI suggests a rate of convergence for h close to 
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N 

Ei[H] 

rate for Ei[H] 

Ei[U] 

rate for Ei [U] 

10 

6.9247 X 10"^ 

- 

1.5110 X 10-2 

- 

20 

3.3950 X 10-2 

1.0283 

2.6966 X 10-3 

2.4863 

40 

1.6779 X 10-2 

1.0168 

5.7095 X 10-1 

2.2397 

80 

8.3767 X 10-3 

1.0022 

1.3474 X 10-1 

2.0831 

160 

4.1860 X 10-3 

1.0008 

3.3114 X 10-3 

2.0247 

320 

2.0924 X 10-3 

1.0004 

8.3275 X 10-® 

1.9915 

640 

1.0460 X 10-3 

1.0002 

2.0570 X 10-® 

2.0173 

Table 8. Spatial 

errors and rates of convergence 

for the exact 

solution with Sh = 

pi and Su = P^ 

using the pi norm. 


N 

Poo[P] 

rate for Poo[P] 

Eoo[U] 

rate for Poo [t^j 

10 

1.3308 X 10-3 

- 

1.4064 X 10-3 

- 

20 

2.6812 X 10-1 

2.3113 

9.9809 X 10-3 

3.8167 

40 

6.2091 X 10-3 

2.1104 

9.1910 X 10-3 

3.4409 

80 

1.5350 X 10-3 

2.0161 

1.0541 X 10-6 

3.1241 

160 

3.8219 X 10-3 

2.0059 

1.2825 X 10-^ 

3.0390 

320 

9.5358 X 10-'^ 

2.0029 

1.6047 X 10-3 

2.9986 

640 

2.3815 X 10-'^ 

2.0014 

1.9751 X 10-3 

3.0223 


Table 9. Spatial errors and rates of convergence for the exact 
solution with Sh = and Su = using the L°° norm. 


3.5 which is in a good agreement with the estimate \\h— ii|| = 0(Aa;^'^-\/ln(l/Aa;)) 
suggested for the standard Galerkin method for Peregrine’s system [3]. Also Table 
[TOl suggests a rate of convergence for h close to 4, which agrees with the estimate 
||m — m| 1 = 0(Aa;'*Y^ln(l/Aa;)) of [4]. 


N 

Eo[H] 

rate for Po[P] Eo[U] 

rate for Eo[U] 

200 

0.5072 X 10-3 

- 

0.3101 X 10-^3 

- 

250 

0.2287 X 10-3 

3.5697 

0.1270 X 10-13 

3.9992 

300 

0.1202 X 10-3 

3.5256 

0.6127 X 10-11 

3.9991 

350 

0.6998 X 10-3 

3.5125 

0.3307 X 10-11 

3.9989 

400 

0.4384 X 10-3 

3.5019 

0.1939 X 10-11 

3.9975 

450 

0.2913 X 10-3 

3.4710 

0.1210 X 10-11 

4.0006 

500 

0.2021 X 10-3 

3.4679 

0.7968 X 10-13 

3.9708 

Table 10. Spatial 

errors and rates of convergence 

for the exact 


solution with finite elements using the norm. 


To study further the accuracy of the method we compute the conservation of 
the energy functional I(t). As initial condition we use a solitary wave of amplitude 
a = 0.2 with a smooth bottom given by the function b{x) = 1 + 0.1 sin(7ra:/2). 
Because the SGN system together with the wall boundary conditions does not 
conserve the energy functional, we consider the interval [—100,100] and maximum 
time T = 50 where the solution remained practically zero at the endpoints of the 
domain. We also used a fairly coarse grid with Ax = 0.1. Due to the dissipative 
properties of the explicit Runge-Kutta method, conservation of the invariant I{t) is 
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N 

Ei[H] 

rate for Ei[H] 

Ei[U] 

rate for Ei [U] 

200 

0.3406 X 10"^ 

- 

0.3920 X lO-"^ 

- 

250 

0.1970 X 10-5 

2.4531 

0.2008 X 10-^ 

2.9982 

300 

0.1265 X 10-5 

2.4294 

0.1162 X 10-^ 

2.9986 

350 

0.8700 X 10-5 

2.4299 

0.7322 X 10-8 

2.9987 

400 

0.6286 X 10-5 

2.4334 

0.4905 X 10-8 

2.9988 

450 

0.4722 X 10-5 

2.4291 

0.3446 X 10-8 

2.9989 

500 

0.3654 X 19-5 

2.4348 

0.2512 X 10-8 

2.9989 

Table 11. Spatial errors and rates of convergence 
solution with hnite elements using the norm. 

for the exact 


N 

E2[H] 

rate for E 2 [H] 

E2[U] 

rate for E 2 [U] 

200 

0.2759 X lO-"" 

- 

0.5098 X 10-^ 

- 

250 

0.2033 X 10-2 

1.3673 

0.3262 X 10-^ 

2.0010 

300 

0.1583 X 10-2 

1.3715 

0.2265 X 10-4 

2.0009 

350 

0.1278 X 10-2 

1.3863 

0.1663 X 10-4 

2.0006 

400 

0.1060 X 10-2 

1.3999 

0.1273 X 10-4 

2.0005 

450 

0.8987 X 10-8 

1.4078 

0.1006 X 10-4 

2.0004 

500 

0.7741 X 10-8 

1.4168 

0.8151 X 10-5 

2.0003 

Table 12. Spatial errors and rates of convergence 
solution with hnite elements using the norm. 

for the exact 


N 

£^00 [77] 

rate for Eoc,[H] 

Eoo[U] 

rate for E^c [U] 

200 

0.5211 X 10-^ 

- 

0.7553 X 10-4'^ 

- 

250 

0.2700 X 10-^ 

2.9460 

0.3111 X 10-45 

3.9751 

300 

0.1577 X 10-^ 

2.9493 

0.1505 X 10-45 

3.9812 

350 

0.1000 X 10-^ 

2.9559 

0.8145 X 10-44 

3.9843 

400 

0.6732 X 10-8 

2.9628 

0.4784 X 10-44 

3.9857 

450 

0.4748 X 10-8 

2.9649 

0.2990 X 10-44 

3.9888 

500 

0.3471 X 10-8 

2.9727 

0.1964 X 10-44 

3.9870 


Table 13. Spatial errors and rates of convergence for the exact 
solution with finite elements using the L°° norm. 


sensitive to the size of the timestep At. For example, when we used Ax = At = 0.1 
with elements in S^, the invariant remained at the value I{t) = 0.314542, when 
we used At = 0.01, then the invariant remained at the value I{t) = 0.31454249795 
conserving the digits shown. 

Although the energy functional I{t) depends on the choice of At, the value of 
the Hamiltonian remains practically constant for any small value of At < 0.01 and 
the error |/(t) —/(0)| is of 0(10“^^). Figured] shows the sensitivity of the invariant 
I{t) to At for fixed Ax = 0.1 for hnite elements. Figure [D shows the error of 
the energy functional as a function of Aa; for different hnite element spaces. We 
observe that the convergence of the computed energy functional is connected to the 
convergence of the numerical method. Although the explicit Runge-Kutta methods 
we used are not conservative (in the sense that they introduce a small amount 
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Figure 1. Sensitivity of energy conservation to the value of At 
for fixed Ax = 0.1 



Ax 


Figure 2. Sensitivity of energy conservation to the value of Ax 
for small values of At 

of numerical dissipation), and their contribution to the error in the conservation 
of energy is larger than the error coming from the spatial discretization of the 
same order, it seems that the actual error is not important and upon choosing 
appropriate small values of At, it can be considered negligible. When we used 
and elements, then the error from the spatial discretisation is larger than the 
errors embedded by the temporal discretisation that are considered unimportant. 

4. Numerical experiments 

In this section we present a series of numerical experiments that serve as bench¬ 
marks to verify the accuracy of the modified Galerkin method, in simulations with 
variable bathymetry and wall boundary conditions. We tested all the numerical 
methods for all the numerical experiments and we report all the important differ¬ 
ences in the numerical results. All the graphs contain the results obtained with P^ 
elements. 


4.1. Shoaling of solitary waves. Shoaling of solitary waves and the nonlinear 
mechanisms behind the shoaling of a solitary wave have been studied theoretically 
and experimentally by various authors, [28l [29l [58l [59], and shoaling is closely 
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related to the runup of solitary wave on a plain beach, m- In |59j after reviewing 
the derivation of Green’s law for the amplitude evolution of shoaling waves for BT 
systems, that is, 

lymaxM ~ , (16) 

several experiments with shoaling of solitary waves are presented. In [58] Green’s 
law rymaxM b{x)~^ for the Shallow Water Wave equations is derived. In [5^ 135] 
four regions of shoaling were determined, where the general Green’s law r^max 
b(x)~°‘ applies for different values of the parameter a: The zone of gradual shoaling 
(a < 1), the zone of rapid shoaling a > 1, the zone of rapid decay a < —1 and the 
zone of gradual decay a > — 1. Green’s law for Boussinesq systems is valid in the 
zone of gradual and rapid shoaling, while Green’s law for the shallow water system 
derived in |58] applies in the zone of rapid shoaling. In this work, we consider 
first the shoaling of solitary waves on a plane beach of (mild) slope 1 : 35. This 
experiment has been proposed by Grilli et. al. [2H1 El] . 

Here we use experimental data taken from |28) to compare with our numerical 
solution. For the numerical simulation we considered the domain shown in Figure 
131 For this experiment we take a uniform mesh on [—100,34] with meshlength Aa: = 
O.I. We also translate the solitary waves so that the crest amplitude is achieved at 
X = —20.1171. Because of the low regularity of the bottom at x = 0, artifacts of 
the numerical method might appear especially in the case of elements. For this 
reason, we ensure the required regularity of the seafloor (required by the model 
equations), by approximating the bottom topography function with a quadratic 
polynomial near x = 0. 

Although the SGN system cannot model accurately the breaking of solitary 
waves, it appears that it models shoaling with higher accuracy than other Boussi¬ 
nesq models [25], even very close to the breaking point. For example when A = 0.2 
the steep wave observed in the laboratory at the location called gauge 9 of [28] is 
approximated by a smooth solution in the SGN system. 

Next, we study the shoaling of solitary waves with normalized amplitude A = 0.1, 
0.15, 0.2 and 0.25. We monitor the numerical solution on the gauges (enumerated as 



X 

Figure 3. Sketch of the domain for the shoaling of solitary waves 
on a plain beach of slope 1 : 35. 

in [28]) with number 0, 1, 3, 5, 7 and 9 located at x = —5.0, 20.96, 22.55, 23.68, 24.68, 25.91 
respectively. The results with elements in the case of the shoaling of the solitary 
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wave with A = 0.2 are shown in Figure IH The results with the other methods are 
almost identical and are not shown here. 




t 


Figure 4. Solution on various wave gauges for the shoaling on a 
plane beach of slope 1 : 35 of a solitary wave with A = 0.2. Circles 
show experimental data, and lines show numerical solutions. 


In these experiments, we also monitored the invariant I(t) for values of t small 
enough that the waves do not interact with the boundaries. In the case of 
elements on a uniform grid with Ax = 0.1, the values of the invariant up to 
time t = 30 are given in the Table [TH For the same experiments, the respec¬ 
tive invariants were conserved to 4 decimal digits in the case of elements with 
Ax = 0.1. Finally, we computed the relative wave height defined as H(x*)/b(x*) 


A 

I(t) 

0.10 

0.104058609813 

0.15 

0.197139475070 

0.20 

0.312548348249 

0.25 

0.449208354485 


Table 14. The conserved invariant I{t) for up to t = 30. 


where H (x*) = maxa;{|C(x, t)|} is the crest amplitude of the wave normalized by the 
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X 


Figure 5. Comparison of the computed relative wave height H/h 
with experimental data of [55] for the shoaling of solitary waves on 
a plane beach of slope 1 : 35, and with Green’s law. Numerical 
solutions, ’o’: Experimental data, Green’s law for Boussinesq 
systems 


local depth h{x*) evaluated at the same point x*. Figure |S| presents the experimen¬ 
tal and numerical data. We observe a very good agreement between the numerical 
and the experimental results. We note that the numerical results obtained by the 
numerical solution of the full water wave problem in [28] fit the experimental data 
(see e.g. Fig. 4 of [28]), as well as the SGN model. 

It is noted that the evolution of the maximum of the solution according to Green’s 
law for Boussinesq systems, [59], represented in Figure[5]by a broken line is initially 
close to the numerical and experimental data, then underestimates the evolution 
of the maximum especially during the shoaling of large amplitude solitary waves. 
This indicates that the solitary waves are in between the zone of rapid and gradual 
shoaling since they are nearly breaking waves. On the other hand Green’s law 
for the nonlinear shallow water system, [58] . which is omitted here, overestimates 
the amplitude of the shoaling solitary wave, likely because of the simplifications 
inherent in the derivation of Green’s law. These results verify the conclusions 
made in [59] about the zone of gradual shoaling where a law that is similar to 
Green’s law r^max ^ is valid. In our case we computed by experimentation 

that 77niax ^ 5“^/^ describes quite well the evolution of the amplitude of a shoaling 
solitary wave on a plane beach of slope I : 35. 

4.2. Reflection of solitary waves on vertical wall. We study also the reflection 
of the solitary waves on a vertical wall. In general, a vertical wall can be modeled 
by assuming that there is no flux through the vertical wall, i.e. the horizontal 
velocity of the fluid on the wall is m = 0. The reflection of a solitary wave at 
a vertical wall is equivalent to the head-on collision of two counter-propagating 
solitary waves of the same shape. Because the head-on collision implies that at the 
center of the symmetric interaction = 0 then one might argue that the symmetric 
head-on collision is governed by different mathematical properties compared to the 
interactions of a solitary wave at a wall. 

In practice, during the interaction of a solitary wave with a wall the derivative rjx 
on the boundary is negligible and so the additional condition seems to be satisfied 
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and the reflection at the wall is apparently equivalent to the head-on collision of 
two counter propagating waves, [3]. (Although the boundary condition rjx = 0 is 
not necessary, its imposition does not require any modification to the numerical 
method or to the finite element spaces). The interaction of solitary waves with a 
vertical wall and the head-on collision of the solitary waves has been studied by 
theoretical and numerical means in |41] . In [41) . an asymptotic solution for the 
maximum runup has been derived, namely, if a is the normalized amplitude Ajh 
of the impinging wave, then the normalized maximum runup .Rmax = R-maxlb is 
approximately 

Rmax ^ 2a -|- ia^ -f ^a^. (17) 

In [IT] , the SGN system with horizontal bottom has been solved numerically using 
a finite difference scheme. Finite difference schemes for the SGN equations have 
certain disadvantages, including the introduction of numerical dispersion or dissi¬ 
pation. In addition it is necessary to use less physical boundary conditions on the 
wall, such as hx = u^x = u — ■^{h^Ux)x = 0i [4:1] . A fourth-order compact finite 
volume scheme was used to solve the SGN with general bathymetry in [Min]- 
The derivation of the boundary conditions proposed in m follow analogous tech¬ 
niques with |48| imposing also more boundary conditions for both h and u that 
approximate the wall boundary conditions. Other useful boundary conditions such 
as absorbing boundary conditions were constructed in m and can be implemented 
equally well with the present numerical model. 

The head-on collision of solitary waves in the full water wave problem has also 
been studied asymptotically, in [56] . Specifically, a similar formula to the one for 
the SGN equations has been derived, and is given by the formula 2a -|- ^a^ -I- |a^. 
Other somewhat less accurate asymptotic approximations have been derived in [46] . 

Recent numerical experiments showed that BT models do not describe accurately 
the reflection of large amplitude solitary waves at a vertical wall. Specifically, in 
[15] solving the full water wave equations using boundary element methods, it was 
shown that large amplitude waves can achieve a higher maximum runup during a 
head-on collision with a vertical wall, than the predicted values of the respective 
asymptotic solutions of [sum] for solitary wave reflection at a vertical wall. These 
results have been verified in [uni]. Additionally, the formation of a residual jet 
was observed during the head-on collision of two large amplitude solitary waves, 
which was responsible for the large values of the maximum runup height. 

To verify the ability of the numerical method to model accurately the reflection 
of solitary waves, we considered the case of a horizontal bottom b{x) = —1 for 
X G [—100, 0] and solitary waves with amplitudes A = 0.075,0.1, 0.15, • • • , 0.7. This 
case and these waves match the experiments which are presented in [18] and which 
serve as benchmarks for the models in [13 [T5]. We also consider a uniform grid on 
the interval [—100,0] with Ax = 0.1 while we translate the solitary waves so that 
their crest amplitude is at a; = —50. Figure[3shows the reflection of solitary waves 
with A = 0.075 and 0.65. 

The reflection at a vertical wall of the extreme case of solitary waves of ampli¬ 
tude A = 0.7 is presented in Figure [T] together with the solution of the head-on 
collision of two symmetric solitary waves. In this figure, we cannot observe any 
differences between the solutions of the reflected wave and of the colliding solitary 
waves (within graphical accuracy). During this interaction a jet-like structure is 
visible, similar to those reported in [H [m HD]- The jet formed in the case of 
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Figure 6. Reflection of solitary waves at a vertical wall located 
at x = 0. 



X 

Figure 7. Reflection of solitary waves of normalized amplitude 
A = 0.7 at a vertical wall located at x = 0, plotted together with 
the head-on collision of two solitary waves of the same amplitude 
to show their equivalence. 


the SGN system is a solution of the mathematical model while the respective jet 
computed in m is not. It is noted that no wave breaking is observed during these 
head-on collisions. Due to the inelastic interaction, the reflected solitary wave is 
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Amplitude 


Figure 8. Maximum runup values during the reflection of solitary 
waves at a vertical wall located at a; = 0. 


followed by a dispersive tail. The magnitude of the dispersive tail depends on the 
amplitude of the colliding solitary wave. As it can be observed in Figure [6] the 
reflection of large solitary waves results in the generation of large dispersive tails. 

In [13] , it was reported that the formation of a residual jet for the full water wave 
equations starts when the normalized amplitude of two incident solitary waves is 
larger than 0.60. In the case of the SGN equations, this jet formation was observed 
in reflections of solitary waves with amplitudes larger than 0.65. A comparison of 
the numerical maximum runup values and the asymptotic formula (El) in Figure 
[5] shows a very close match between the numerical solutions and the asymptotic 
results. In this figure, the numerical results of both the piecewise linear and the 
cubic spline finite elements are presented but no differences can be observed within 
graphical accuracy. Moreover, the numerical results of [18] are compared with 
the numerical results obtained with the modified Galerkin method and we verify 
the difference in the head-on collision processes between the SGN and the Euler 
equations, [Ullla]. This can be explained by noting that the jet cannot be modeled 
by the SGN equations in the present form, since this artifact cannot be described 
by a smooth function but only by a parametric curve, m- 

Remark 4.1. In some cases it is considered that the reflection of a wave-train of 
more than one pulse can result in extreme runup values. Extreme wave runup on 
a vertical wall has been studied using periodic-boundary conditions and the head-on 
collision of wave-trains propagating in different directions in [60] [12]. We repeated 
several of these experiments using half of the domain (reguired by the periodic code 
and with the wall boundary conditions described in this paper). The results were 
identical (except for the accuracy) to those obtained in [50] Ej while our code re¬ 
mained stable during the strong interactions of the waves with the wall. This verifies 
that the wall boundary conditions and the new numerical method can describe ac¬ 
curately strong interactions and extreme wave runup on a wall. As an indication 
of the results obtained we report here only the case of the reflection of a wave-train 
with three pulses = 3, Aq = 125c? in the notation of [60 ] ). In this case the 
normalized maximum runup observed was i?max/ao ~ 5.545. 









A MODIFIED GALERKIN METHOD EOR THE SGN SYSTEM 


23 


4.3. Reflection of shoaling waves. In our final set of numerical experiments we 
consider the benchmark described in [ 6 ll |20] . Specifically, we consider the compu¬ 
tational domain [—100, 20] (with Ax = 0.1) and a bottom topography consisting of 
a horizontal seafloor in [—100,0] and a plane beach of slope 1 : 50 for a; £ [0, 20]. A 
sketch of the computational domain is presented in Figure |9l The vertical wall is 
located at x = 20 . 

In this set of experiments we study the reflection of solitary waves at a vertical 
wall after the waves have climbed up the sloping beach. We consider two cases, one 
with a solitary wave of amplitude A = 0.07 and another with A = 0.12. The initial 
conditions have been translated so that the maximum of the crest is at x = —30. 
Figures dO] and [H] show experimental data recorded at three locations, by wave 
gauges gi, 52 , and 53 , placed at x = 0, 16.25 and 17.75. The comparison of numer¬ 
ical solutions with experimental data shows a better match from the present SGN 
model than is obtained using any weakly nonlinear and weakly dispersive Boussi- 
nesq system, |25l 161) . These numerical experiments verify the ability of the present 
numerical scheme to approximate with high accuracy and model fully nonlinear 
and weakly dispersive waves with wall boundary conditions at the endpoints of the 
computational domain. 



Figure 9. Sketch of the numerical experiment for the reflection 
at a vertical wall located at x = 20 of a shoaling wave over a plain 
beach of slope 1 : 50. 

In order to study the stability of the modified Galerkin method in more demand¬ 
ing situations, we consider the propagation of a solitary wave over a composite 
beach simulating the geometrical dimensions of Revere Beach, and its reflection 
by vertical wall. These experiments were conducted at the Goastal Engineering 
Laboratory of the U.S. Army Corps of Engineers, Vicksburg, Mississippi facility, 
[33] and serve as benchmarks for the reflection of nonbreaking, nearly breaking and 
breaking solitary waves by vertical wall. The composite beach consists of three 
piecewise linear segments while the bathymetry is constant away from the beach 
and equal to bo = 0.218 m. The bathymetry can be realized by the function: 

11.77 < X < 15.04 
15.04 < X < 19.4 
19.4 < X < 22.33 ■ 

22.33 < X < 23.23 


b{x) = 


-0.218, 

1/53 X- 0.5018, 
1/150 X- 0.2650, 
1/13 X- 1.8340, 
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Figure 10. Reflection at a vertical wall located at a: = 20 of a 
shoaling wave over a plain beach of slope 1 : 50. Initial solitary 
wave amplitude A = 0.07. 


The sketch of the domain is presented in Figure [121 Here, we considered 
three solitary waves with normalized amplitudes A/bo = 0.05, 0.3 and 0.7. We 
monitored the water depth at x locations that correspond to gauges 5, 7 and 9 
in [33]. In this experiment, we used Ax = 0.1 in [—11.77,23.23]. The normalized 
maximum runup computed at the location of the vertical wall in the first case 
was R/bo = 0.122 which is very close to the experimental value measured in the 
laboratory R/bo = 0.13. Moreover, the computed solution is very close to the 
experimental data at the gauges and we don’t present these results here. It is 
noted that because a very small relative amplitude wave is involved in this case, 
its reflection can be modelled quite accurately, even by nondispersive models. The 
other two cases involve a nearly breaking and a breaking wave. Due to the steepness 
of the wave in the third case a wave breaking mechanism should be considered in 
order to approximate the solution in a stable manner. In the second experiment, 
the solitary wave is a nearly breaking wave. Although in this case the wave becomes 
very steep during shoaling, the numerical maximum runup computed was 0.46 to, 
which is again very close to the experimentally recorded runup value 0.45 to. The 
solution at the wave gauges is presented in Figure [T31 The results can be improved 
by considering wave breaking mechanisms |541 [7l [64] and Green-Naghdi models 
with improved dispersion characteristics such as those proposed in [ini[Il]|6][5][37]. 
It is noted that in both cases only linear terms with second order derivatives are 
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Figure 11. Reflection at a vertical wall located at a; = 20 of a 
shoaling wave over a plain beach of slope 1 : 50. Initial solitary 
wave amplitude A = 0.12. 



X 


Figure 12. Sketch of the domain for the reflection of a solitary 
wave over a composite beach. 


included locally around the regions where the waves become steep. These terms can 
be incorporated easily by any numerical method while the possible stiffness induced 
by the new terms to the problem can be handled by taking smaller time-steps or 
by using a time integration method appropriate for stiff problems, [7]. 
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Figure 13. Reflection of a solitary wave over a composite beach. 
Vertical wall located at 23.23. Initial solitary wave amplitude 
A/6o = 0.3. 


5. Summary and Conclusions 

We present a fully discrete numerical scheme for the Serre-Green-Naghdi (SGN) 
system with wall boundary conditions. Semidiscretization of the model equations is 
based on a modification of the standard Galerkin / finite-element method that al¬ 
lows solutions in function spaces of low regularity. The time discretization is based 
on a fourth-order, four-stage, explicit Runge-Kutta method. A detailed computa¬ 
tional study of the convergence properties of the numerical scheme shows that the 
method converges with similar convergence rates to those of Peregrine’s system. 
Some of the advantages of the new numerical method are: 

• the high accuracy and the very good conservation properties; 

• the sparsity of the resulting linear systems; 

• the low complexity of the algorithm due to the use of low-order finite ele¬ 
ment spaces; and 

• its potential to be extended to the two-dimensional model equations. 

In addition, we perform a series of very accurate numerical experiments to verify 
the efficacy of the numerical scheme in studies of shoaling and reflecting solitary 
waves. The numerical solutions are compared with available experimental data 
and theoretical approximations, whenever possible. The numerical model appears 
to be efficient and the match between numerical results, experimental data, and 
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theoretical approximations is very satisfactory and shows better performance than 
other shallow water wave systems when there is no wave breaking. Wave breaking 
can be treated by following heuristic methodologies, such as discarding dispersive 
terms or adding new dissipative terms. 
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